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INTRODUCTION: 

A  tracked  free-hand  ultrasound  system  is  ideal  for  guiding  many  radiotherapy  procedures,  as  it  can  be  performed  in 
the  treatment  room.  In  partial  breast  irradiation,  the  lumpectomy  cavity  should  be  localized  during  the  treatment 
course.  While  many  structures  in  the  breast  look  similar  to  the  lumpectomy  cavity  in  ultrasonography,  the  scar  tissue 
around  the  cavity  is  hard  and  can  be  visualized  by  ultrasound  elastography.  To  be  clinically  successful,  the 
elastography  method  needs  to  be  robust  to  the  sources  of  decorrelation  between  ultrasound  images,  specifically  fluid 
motions  inside  the  cavity,  change  of  the  appearance  of  speckles  caused  by  compression  and  physiologic  motions, 
and  out-of-plane  motion  of  the  probe.  In  this  paper,  we  present  a  novel  elastography  technique  that  is  based  on 
analytic  minimization  of  a  regularized  cost  function.  The  cost  function  incorporates  similarity  of  RF  data  intensity 
and  displacement  continuity,  making  the  method  robust  to  decorrelation  noise  present  throughout  the  image.  We 
exploit  techniques  from  robust  statistics  to  make  the  method  resistant  to  large  decorrelations  caused  by  sources  such 
as  fluid  motion.  The  analytic  displacement  estimation  works  in  real-time.  The  tracked  data,  used  for  targeting  the 
irradiation,  is  also  exploited  for  discarding  frames  with  excessive  out-of-plane  motion.  Simulation,  phantom  and 
patient  results  are  presented. 


BODY: 

We  are  very  excited  to  report  the  novel  achievements  of  this  research  effort.  The  detailed  Statement  of  Work  tasks 
and  a  description  below  each  task  is  provided  next.  The  SOW  tasks  are  in  blue.  A  description  of  the  research  effort 
follows  each  SOW  item  in  black. 


1.  Obtain  ultrasound  (US),  US  elastography  (USE)  and  CT  scans  before  the  start  of  the  radiotherapy. 

a.  Develop  a  tracked  US  system  for  data  collection  (month  1).  We  will  use  Polaris  optical  tracker, 
already  available  in  the  ERC-CISST  for  tracking.  Polaris  markers  will  be  attached  to  the  US  probe 
and  the  US  probe  will  be  calibrated  with  respect  to  the  trackers. 

Polaris  optical  tracker  gives  very  accurate  displacement  measurements.  However,  they  require  line  of 
sight  (as  they  are  optical)  and  also  they  take  longer  time  to  set  up  in  a  CT  imaging  room.  Because  of  the 
time  constraints  (we  had  a  short  amount  of  time  to  set  up  the  tracking  device  in  the  CT  room  before  the 
patient  shows  up),  we  instead  used  magnetic  trackers.  Magnetic  trackers  are  not  as  accurate  as  Polaris 
and  suffer  from  distortion  of  magnetic  field  due  to  presence  of  metals.  However,  they  do  not  require 
line  of  sight  and  also  easier  to  set  up  and  also  are  inexpensive. 


b.  Collect  US  data  from  patient  before  the  PBI  treatment  at  the  same  time  that  CT  is  collected 
(months  2-14). 
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We  have  collected  data  from  more  than  10  patients  so  far  and  continue  to  acquire  more  data.  All  data  is 
acquired  in  the  CT  room  while  the  patient  is  lying  down  on  the  CT  bed.  The  patient  does  not  move 
after  ultrasound  data  is  acquired  until  her  CT  scan  is  finished.  We  have  acquired  both  3D  ultrasound 
data  (freehand  data,  no  3D  probe  is  used)  and  palpation  data  (for  elastography).  More  details  of  the 
system  are  in  the  attached  paper  (MICCAI  conference,  34%  acceptance  rate,  MedLine  and  ISI  indexed 
and  listed).  I  received  travel  award  from  MICCAI  (only  7  students  received  this  award)  to  present  this 
work. 


c.  Register  US  to  the  CT  (months  2-14). 

The  registration  of  ultrasound  to  CT  is  underway.  We  have  so  far  looked  at  the  data  which  is  available 
from  the  tracking  device  (the  magnetic  tracker)  to  aid  the  registration.  Development  of  a  non-rigid 
registration  technique  that  registers  CT  and  ultrasound  using  mutual  information  technique  (an  image 
based  registration  technique)  is  underway.  Registration  of  ultrasound  to  elastography  is  pretty  straight 
forward  and  has  been  completed. 


d.  Compare  different  combination  of  US,  USE  and  CT  for  delineation  (months  3-18). 

We  have  done  preliminary  studies  on  comparing  the  performance  of  the  operator  with  different  imaging 
modalities.  While  it  is  intuitive  that  adding  ultrasound  and  elastography  enhance  the  delineation  of  the 
lumpectomy  bed,  we  have  not  achieved  a  conclusive  result  yet. 


e.  Optimize  USE  code  for  using  the  human  data  (months  3-24).  I  have  developed  a  novel  USE 
technology  and  software  under  a  Breast  Cancer  Research  Foundation  research.  The  work  thus  far 
has  been  extremely  promising;  a  manuscript  has  already  been  accepted  for  publication  in  the  IEEE 
Trans.  Med.  Imag.  I  will  be  further  improving  my  USE  implementation  to  enhance  visualization  of 
the  lumpectomy  bed  and  ductal  tissue. 

Extremely  promising  results  have  been  obtained  in  this  area.  We  have  developed  new  elastography 
techniques  that  generate  superb  images  from  the  human  data.  The  attached  MICCAI  paper  contains 
some  of  the  details  of  the  method.  Another  journal  paper  is  submitted  to  the  IEEE  TMI  (the  journal 
with  the  highest  impact  factor  in  medical  imaging).  Since  the  paper  is  still  in  the  submission  stage,  we 
cannot  attach  it  to  this  report.  It’s  reference  is: 

Rivaz,  H.,  Boctor,  E.,  Choti,  M.,  Hager,  G.,  “Real-time  Regularized  Ultrasound  Elastography”,  IEEE 
Trans.  Medical  Imaging  (submitted) 


f.  Performance  optimization  and  refinements  of  subsystems  (months  18-36). 

The  ultrasound  elastography  subsystem  is  almost  finalized.  We  have  made  the  code  also  available  to 
the  public,  so  that  groups  who  work  in  different  applications  of  ultrasound  elastography  can  exploit  our 
novel  and  high-quality  elastography  estimation  technique.  We  could  confidently  say  that  the  method 
can  be  implemented  commercially  with  small  modofications.  Work  on  registration  subsystem  is  now 
underway. 
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2.  Obtain  tracked  US  and  USE  scans  of  the  patients  weekly  during  the  radiotherapy  (months  2-14).  This 
data  will  be  acquired  from  the  same  patients  from  whom  US  and  USE  data  were  acquired  in  step  l.b. 

As  mentioned  in  the  lb  Section,  elastography  and  ultrasound  images  are  both  acquired  from  the  patients. 


3.  Develop  US  to  US  registration  (month  25-27).  This  task  is  partly  solved  in  task  l.c.  (US  to  CT 
registration).  As  mentioned,  US  to  US  registration  is  one  of  the  steps  required  in  US  to  CT  registration.  I 
will  optimize  this  step  to  work  with  breast  images,  as  opposed  to  simulation  images  in  task  l.c. 

Preliminary  work  on  the  theoretical  aspects  of  this  task  has  started. 


4.  Register  the  US  images  to  the  US  images  obtained  before  the  therapy  (months  27-36).  US  images 
obtained  in  task  l.b  have  the  lumpectomy  bed  delineated.  US  images  obtained  in  task  2  are  acquired  from 
the  same  patients,  and  therefore  the  delineated  volume  can  be  projected  into  them  (without  having  to 
delineate  the  lumpectomy  bed  again).  Having  the  position  of  the  lumpectomy  bed  in  the  tracked  US  images 
allows  updating  linear  accelerator  (linac)  guidance.  This  is  important  since  the  position  of  the  lumpectomy 
bed  changes  during  the  long  therapy  session. 

Work  has  not  started  yet. 

KEY  RESEARCH  ACCOMPLISHMENTS: 

•  Development  of  a  novel,  real-time,  robust  and  high  quality  elastography  techniques 

•  Development  of  a  tracked  ultrasound  system  using  magnetic  trackers 

•  Study  of  ultrasound  images  of  lumpectomy  cavity 

•  Study  of  elastography  images  of  lumpectomy  cavity 

•  Data  acquisition  from  more  than  10  patients 

REPORTABLE  OUTCOMES: 

1.  Rivaz,  H.,  Foroughi,  P.,  Fleming,  I.,  Zellars,  R.,  Boctor,  E.,  Hager,  G.,  “Tracked 

Regularized  Ultrasound  Elastography  for  Targeting  Breast  Radiotherapy”,  Medical 

Image  Computing  and  Computer  Assisted  Intervention,  MICCAI,  London,  UK,  Sept. 

2009,  pp  507-515.  [Acceptance  rate:  34%]  [Awarded  MICCAI  Travel  Grant] 

2.  Rivaz,  H.,  Boctor,  E.,  Choti,  M.,  Hager,  G.,  “Regularized  Ultrasound  Elastography”,  IEEE  Trans.  Medical 

Imaging  (highest  impact  factor  in  medical  imaging  engineering  journals).  Submitted. 

3.  Rivaz,  H.,  Fleming,  I.,  Assumpcao,  L.,  Fichtinger,  G.,  Hamper,  U.,  Choti,  M.,  Hager,  G.,  Boctor,  E.,  “Ablation 

Monitoring  with  Elastography:  2D  In-vivo  and  3D  Ex-vivo  Studies”,  Medical  Image  Computing  and  Computer 

Assisted  Intervention,  MICCAI,  New  York,  NY,  Sept.  2008,  pp  458-466  .  [Acceptance  rate:  35%] 


CONCLUSION: 

Data  acquisition  from  lumpectomy  patients  went  smoothly.  We  acquired  the  ultrasound  data  few  minutes  before  the 
CT  scan  was  acquired.  This  means  that  we  did  not  alter  the  current  medical  work-flow.  Ultrasound  data  acquisition 
took  only  few  minutes,  minimizing  the  cost  and  patient  discomfort.  Magnetic  tracking  provided  ease  of  use  and 
enough  reliability  for  registration  of  ultrasound  to  CT  images  and  for  reconstructing  3D  volumes  from  2D  data.  The 
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novel  ultrasound  elastography  technique  has  many  advantages  compared  to  the  previous  methods.  We  chose  the 
novel  application  of  the  lumpectomy  cavity  localization  as  the  hard  scar  tissue  around  the  lumpectomy  is  relatively 
thin  and  demands  a  high  resolution  elastography  method.  Also,  incoherent  fluid  motions  in  the  cavity  causes  large 
decorrelations,  requiring  a  robust  method.  The  merit  of  adding  ultrasound  and  elastography  to  delineation  of  the 
lumpectomy  bed  is  not  clear  to  us  yet.  We  are  performing  more  studies  to  verify  the  importance  of  ultrasound  and 
ultrasound  elastography. 

1.  Rivaz,  H.,  Stolka,  P.,  Hager  G.,  Boctor,  E.  “Robust  sensorless  freehand  3D  ultrasound”,  IEEE  Trans.  Medical 
Imaging  (in  preparation) 

2.  Rivaz,  H.,  Boctor,  E.,  Choti,  M.,  Hager,  G.,  “Regularized  Ultrasound  Elastography”,  IEEE  Trans.  Medical 
Imaging  (submitted) 

3.  Rivaz,  H.,  Shinagawa,  Y.,  Liang,  J.  “Physical  Priors  in  Electronic  Colon  Cleansing  for  Virtual  Colonoscopy”, 

Med.  Phys.  (submitted) 

4.  Rivaz,  H.,  Boctor,  E.,  Eoroughi,  P.,  Zellars,  R.,  Eichtinger,  G.,  Hager,  G.,  “Ultrasound  Elastography:  a  Dynamic 
Programming  Approach”,  IEEE  Trans.  Medical  Imaging  Oct.  2008,  vol.  27  pp  1373-1377 

5.  Rivaz,  H.,  Rohling,  R.,  “An  Active  Dynamic  vibration  Absorber  for  a  Hand-Held  Vibro-Elastography  Probe,” 
ASME  Trans.  Vibration  &  Acoustics,  Eeb.  2007,  vol.  129,  pp  101-112  [Top  10  most  downloaded  articles  from  this 
journal  in  March  2008] 

6.  Rivaz,  H.,  Kang,  H.,  Stolka,  P.,  Boctor,  E.  “Novel  reconstruction  and  feature  exploitation  techniques  for 
sensorless  freehand  3D  ultrasound”,  SPIE  Med.  Imag.,  2010  (accepted) 

7.  Rivaz,  H.,  van  Vledder,  M.,  Choti  M.,  Hager,  G.,  Boctor,  E.  “Liver  ablation  guidance:  discriminating  ablation 
tumor  from  the  cancer  tumor  with  ultrasound  elastography”,  SPIE  Med.  Imag.,  2010  (accepted) 

8.  Rivaz,  H.,  Eoroughi,  P.,  Eleming,  I.,  Zellars,  R.,  Boctor,  E.,  Hager,  G.,  “Tracked  Regularized  Ultrasound 
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9.  Rivaz,  H.,  Liang,  J.,  Shinagawa  Y.,  “Electronic  Colon  Cleansing  of  the  Unprepared  Colon”,  SPIE  Med.  Imag., 
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Tracked  Regularized  Ultrasound  Elastography 
for  Targeting  Breast  Radiotherapy 


Hassan  Rivaz,  Pezhman  Foroughi,  loana  Fleming,  Richard  Zellars, 
Emad  Boctor,  and  Gregory  Hager 

Johns  Hopkins  University,  Baltimore,  MD,  USA 


Abstract.  Tracked  nltrasound  elastography  can  be  used  for  guidance  in 
partial  breast  radiotherapy  by  visualizing  the  hard  scar  tissue  around  the 
lumpectomy  cavity.  For  clinical  success,  the  elastography  method  needs 
to  be  robust  to  the  sources  of  decorrelation  between  ultrasound  images, 
specifically  fluid  motions  inside  the  cavity,  change  of  the  appearance  of 
speckles  caused  by  compression  or  physiologic  motions,  and  out-of-plane 
motion  of  the  probe.  In  this  paper,  we  present  a  novel  elastography  tech¬ 
nique  that  is  based  on  analytic  minimization  of  a  regularized  cost  func¬ 
tion.  The  cost  function  incorporates  similarity  of  RF  data  intensity  and 
displacement  continuity,  making  the  method  robust  to  small  decorre¬ 
lations  present  throughout  the  image.  We  also  exploit  techniques  from 
robust  statistics  to  make  the  method  resistant  to  large  decorrelations 
caused  by  sources  such  as  fluid  motion.  The  analytic  displacement  esti¬ 
mation  works  in  real-time.  Moreover,  the  tracked  data,  used  for  targeting 
the  radiotherapy,  is  exploited  for  discarding  frames  with  excessive  out- 
of-plane  motion.  Simulation,  phantom  and  patient  results  are  presented. 


1  Introduction 

Breast  irradiation  after  lumpectomy  significantly  reduces  the  risk  of  cancer  re¬ 
currence.  There  is  growing  evidence  suggesting  that  irradiation  of  only  the  in¬ 
volved  area  of  the  breast,  partial  breast  irradiation  (FBI),  is  as  effective  as  whole 
breast  irradiation  [1].  Benehts  of  PBI  include  signihcantly  shortened  treatment 
time  and  fewer  side  effects  as  less  tissue  is  treated.  However,  these  benehts  cannot 
be  realized  without  localization  of  the  lumpectomy  cavity.  Tracked  ultrasound 
elastography  can  be  used  for  localizing  the  lumpectomy  cavity  in  the  treatment 
room,  minimizing  tissue  motion  from  planning  to  treatment. 

This  paper  is  focused  on  freehand  palpation  elastography,  which  involves  es¬ 
timating  the  displacement  held  of  the  tissue  undergoing  slow  compression.  Most 
elastography  techniques  estimate  the  displacement  held  using  local  cross  corre¬ 
lation  analysis  of  echoes  [2,3,4].  These  methods  are  very  sensitive  and  accurate 
for  calculating  small  displacements.  However,  elastography  is  subject  to  speckle 
decorrelation  caused  by  various  sources  such  as  motion  of  subresolution  scatter- 
ers,  out-of-plane  motion,  high  compression  and  complex  fluid  motions. 

The  prior  of  tissue  deformation  continuity  can  be  used  to  make  elastography 
more  robust  to  signal  decorrelation.  Previous  work  on  regularized  elastography 
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is  computationally  expensive  [5,6].  Dynamic  programming  (DP)  can  be  used 
to  speed  the  optimization  procedure  [7],  but  it  only  gives  integer  displacements. 
Subpixel  displacement  estimation  is  possible  [7] ,  but  it  is  computationally  expen¬ 
sive  if  a  fine  subpixel  level  is  desired.  In  addition,  a  fixed  regularization  weight  is 
applied  throughout  the  image.  However,  while  two  ultrasound  images  may  cor¬ 
relate  well  in  most  parts,  they  can  have  small  correlation  in  specific  parts.  Four 
examples  of  low  correlation  are:  (1)  correlation  decreases  with  depth  mainly  due 
to  a  decrease  in  the  ultrasonic  signal  to  noise  ratio,  (2)  correlation  is  low  close 
to  arteries  due  to  complex  motion  and  inside  vessels  due  to  blood  motion,  (3) 
correlation  is  extremely  low  in  lesions  that  contain  liquid  due  to  the  incoherent 
fluid  motion  [8,3],  and  (4)  out-of-plane  motion  of  movable  structures  within  the 
image  [8]  causes  low  local  correlation.  To  prevent  such  regions  from  introduc¬ 
ing  errors  in  the  displacement  estimation  one  should  use  large  weights  for  the 
regularization  term,  resulting  in  over-smoothing. 

Freehand  palpation  elastography  provides  ease-of-use  and  requires  minimum 
additional  cost.  However,  out-of-plane  motion  cannot  be  avoided  in  freehand  pal¬ 
pation,  which  reduces  the  quality  of  any  elastography  method.  Assisted  freehand 
elastography  [9]  significantly  reduces  the  out-of-plane  motion  but  it  requires  ad¬ 
dition  of  a  device  to  the  probe.  Quality  metrics  such  as  persistence  in  strain 
images  have  also  been  developed  to  address  this  problem  [10].  To  measure  the 
persistence,  elastography  is  performed  on  two  pairs  of  images  and  the  resulting 
strain  images  are  correlated.  This  method  requires  strain  images  for  calculating 
the  quality  metric.  Therefore,  trying  all  the  combinations  in  a  series  of  frames 
to  find  the  best  pair  for  elastography  will  be  computationally  expensive. 

In  this  paper,  we  present  a  novel  elastography  method  based  on  analytic 
minimization  (AM)  of  a  cost  function  that  incorporates  similarity  of  echo  ampli¬ 
tudes  and  displacement  continuity.  We  introduce  a  novel  regularization  term  and 
demonstrate  that  it  minimizes  displacement  underestimation  caused  by  smooth¬ 
ness  constraint.  We  also  introduce  the  use  of  robust  statistics  implemented  via 
iterated  reweighted  least  squares  (IRLS)  to  treat  uncorrelated  ultrasound  data 
as  outliers.  And  finally,  we  use  the  tracking  information  to  select  the  best  pairs 
of  frames  for  elastography.  Simulation,  phantom  and  patient  experiments  are 
presented  for  validation. 

2  Regularized  Displacement  Estimation 

Dynamic  Programming  (DP).  DP  is  a  discrete  efficient  optimization  tech¬ 
nique  for  causal  systems.  In  DP  elastography  [7],  a  cost  function  is  defined  as 

C{i,  di)  =  min{C'(i  -  l,di-i)  -H  aaR{di,di-i)} +  \h{i)  -  hii  +  di)\  ,  i  =  2-  ■  -  m 

di-i 

(1) 

where  di  is  the  displacement  of  sample  z,  R{di^di-i)  —  {di  —  di-i)  is  an  axial 
regularization  term  (axial,  lateral  and  out-of-plane  directions  are  respectively 
2:,  X  and  y  in  Figure  2  (a)),  aa  is  a  weight  for  the  regularization,  R  and  R 
are  corresponding  RF-lines  of  before  and  after  deformation  and  m  is  the  length 
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of  RF-lines.  The  cost  function  is  minimized  at  z  =  m  and  the  di  values  that 
have  minimized  the  cost  function  are  traced  back  io  i  —  1,  giving  the  di  for  all 
samples.  We  have  implemented  a  2D  DP  algorithm  similar  to  [7]  to  generate 
integer  displacements  as  a  starting  point  for  the  next  step  of  our  algorithm. 

Analytic  Minimization  (AM).  We  now  propose  a  method  that  analytically 
minimizes  a  regularized  cost  function  and  gives  the  rehned  displacement  held. 
Only  axial  displacements  will  be  rehned  for  strain  calculation. 

Having  the  integer  displacements  di  from  DP,  it  is  desired  to  hnd  Adi  values 
such  that  di  +  Adi  gives  the  value  of  the  displacement  at  the  sample  i  ior  i  = 
1  ■  ■  -  m.  Such  Adi  values  will  minimize  the  following  regularized  cost  function 

C  {Adi,  ■  ■  •  5  Adm)  =  [hi'i)  —  d2{'i  +  di  +  Adi)]^  + 

aa{di  +  Adi  -  di-i  -  Adi-if  +  ai{di  +  Adi  -  df  -  Adff  (2) 

where  superscript  p.  refers  to  the  previous  RF-line  (adjacent  RF-line  in  the 
lateral  direction)  and  ai  is  a  weight  for  lateral  regularization.  Substituting  l2{i  + 
di  +  Adi)  with  its  hrst  order  Taylor  expansion  approximation  around  di,  we  have 

O  {Adi,  ■  ■  ■ ,  Adm)  =  S^i  [Ii{i)  -l2{i  +  d,)  -  I'S  +  d,)Adi)f  + 

aa{di  +  Adi  -  di-i  -  Adi-i)"^  +  ai{di  +  Adi  -  df  -  Adf)'^  (3) 

where  is  the  derivative  of  the  l2-  The  optimal  Adi  values  occur  when  the 
partial  derivative  of  C  w.r.t.  Adi  is  zero.  Setting  =  0  we  have 

1  -1  0  ■■■O' 
-1  2  -1  ■■■  0 

0  ■■■  0  -11 
(4) 

rp  rp 

where  I2  =  diag{l2{l  +  di)  ■  ■  ■  /^(m  +  d^n)),  Ad  =  [Z\di  ■  ■  ■  Adm]  ,  e  =  [ei  ■  ■  ■  Cm]  , 
ei  =  Ii{i)  —  l2{i  +  di),  d  =  [di  ■  ■  ■  d^]^,  d^'^-  =  +  Ad^  is  the  vector  of  total 

displacement  of  the  previous  line  and  I  is  the  identity  matrix.  I2,  D  and  I  are 
matrices  of  size  m  x  m  and  Ad,  r,  d  and  d^'^'  are  vectors  of  size  m. 

Biasing  the  Regularization.  The  regularization  term  aa{di  +  Adi  —  di-i  — 
Adi-i)^  penalizes  the  difference  between  di+Adi  and  di-i+Adi-i,  and  therefore 
can  result  in  underestimation  of  the  displacement  held.  Such  underestimation 
can  be  prevented  by  biasing  the  regularization  by  e  to  aa{di  +  Adi  —  di-i  — 
Adi-i  —  e)^,  where  e  =  {dm  —  di)/(m  —  1)  is  the  average  displacement  difference 
between  samples  i  and  i  —  1.  An  accurate  enough  estimate  of  dm  —  di  is  known 
from  the  previous  line.  With  the  bias  term,  the  R.H.S.  of  Equation  4  becomes 
I^e  —  (ctaD  +  CK/I)d  +  Q!;d*-P'  +  b  where  the  bias  term  is  b  =  aa[—€  0  ■  ■  ■  0  e] 
and  all  other  terms  are  as  before.  Interestingly,  except  for  the  hrst  and  the  last 
equation  in  this  system,  all  other  m  —  2  equations  are  same  as  Equation  4. 

Equation  4  can  be  solved  for  Ad  in  4m  operations  since  the  coefficient  matrix 
I2  +  c^aD  +  CT/I  is  tridiagonal.  Utilizing  its  symmetry,  the  number  of  operations 


(I^"  +  a„D  +  ad)  Ad  =  I^e  -  (a^D  +  ad)d  +  a/d^'P',  D  = 
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can  be  reduced  to  2m.  The  number  of  operations  required  for  solving  a  system 
with  a  full  coefficient  matrix  is  more  than  m^/3,  signihcantly  more  than  2m. 

Making  Tracking  Resistant  to  Ontliers.  Even  with  pure  axial  compression, 
some  regions  of  the  image  may  move  out  of  the  imaging  plane  and  increase 
the  decorrelation.  In  such  parts  the  conhdence  of  the  data  term  is  less  and 
therefore  the  weight  of  the  regularization  term  should  be  increased.  The  parts 
of  the  image  with  low  correlation  can  be  regarded  as  outliers  and  therefore 
a  robust  estimation  technique  can  limit  their  effect.  Before  deriving  a  robust 
estimator  for  Ad,  we  rewrite  Equation  3  as  C{Ad)  =  U'^]^p{ri)  +i?(Z\d)  where 
Ti  =  Ii{i)  —  l2{i  +  di)  —  l2{i  +  di)Adi,  p{ri)  =  rf  and  R  is  the  regularization  term. 
The  M-estimate  of  Ad  is  Ad  =  axgm.m.Ad{R^iP{Ti)  +i?(/Ad)}  where  p{u)  is 
a  robust  loss  function  [11].  The  minimization  is  solved  by  setting  =  0: 


dr  dR{Ad) 
^^"’dAdi  dAd^ 


0 


(5) 


A  common  next  step  [11]  is  to  introduce  a  weight  function  w,  where  w{ri).ri  — 
p'iri).  This  leads  to  a  process  known  as  “iteratively  reweighted  least  squares” 
(IRLS),  which  alternates  steps  of  calculating  weights  w{ri)  for  =  1  ■  ■  -  m  using 
the  current  estimate  of  Z\d  and  solving  Equation  5  to  estimate  a  new  Z\d  with 
the  weights  hxed.  Among  many  proposed  shapes  for  w{-),  we  use  [11] 


w{ri) 


1 


n\  <  T 
ri\  >  T 


(6) 


where  T  is  a  threshold  that  can  be  tuned.  A  small  T  will  treat  many  samples  as 
outliers.  With  the  addition  of  the  weight  function,  Equation  5  becomes 

(wl2^  +  ctD  +  CK2i)  Ad  =  wl^e  —  (ctiD  +  Q:2i)d  +  CK2d*'^'  +  b  (7) 


where  w  =  diag{w{ri)  ■  ■  ■w{rm))-  All  of  the  results  presented  in  this  work  are 
obtained  with  one  iteration  of  the  above  equation  unless  otherwise  specihed. 
Current  implementation  of  the  AM  algorithm  with  the  IRLS  takes  0.015s  to 
generate  a  dense  displacement  held  of  size  1300  x  60  on  a  3.4GHz  P4  CPU(not 
including  the  DP  run  time).  The  computation  time  increases  linearly  with  the 
size  of  images. 

Frame  Selection.  The  ultrasound  probe  is  tracked  in  navigation/guidance 
systems  to  provide  spatial  information,  to  generate  freehand  3D  ultrasound, 
or  to  facilitate  multi-modality  registration.  Through  a  calibration  process,  the 
6DOF  motion  of  the  probe  in  the  sensor  coordinate  system  is  transformed  into 
image  coordinate  system  [12].  The  mean  of  the  absolute  motion  value  of  all  pixels 
in  3D,  {|n3;|),  (|ny|)  and  (|n^|),  can  be  analytically  related  to  the  6DOF  sensor 
readings  using  straightforward  and  efhcient  geometric  computations.  For  frame 
i  and  j  to  be  selected  from  a  sequence  of  frames  for  elastography. 


{\Vz\)  -Vz,opt\ 
{\Vz\)+C 


Qi,i  =  K  (|na:|)^  -t-  ky  {\Vy\f  -t-  /c. 


3 


(8) 
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should  be  minimized  where  kx-,  ky,  and  kz  are  weights  for  lateral,  out-of-plane 
and  axial  displacements  and  Vz,opt  is  the  optimum  axial  motion.  Please  refer  to 
[12]  for  a  rationale  of  the  shape  this  function.  Note  that  the  selected  pairs  are 
not  necessarily  consecutive  frames.  The  parameters,  kx,  ky,  kz,  Vz,opt  and  c  are 
manually  tuned  to  1,  2,  1,  0.7  and  1  for  the  AM  elastography  method. 

3  Simulation,  Phantom  and  Patient  Results 

Simulation  Results.  RF  ultrasound  data  of  two  phantoms  are  simulated  using 
Field  II  [13].  The  hrst  phantom  is  50  x  10  x  55mm  and  the  second  one  is  36  x 
10  X  25mm.  They  are  both  made  of  homogeneous  and  isotropic  material:  the 
hrst  one  is  uniform  and  the  second  one  contains  a  circular  hole  hlled  with  water, 
simulating  a  blood  vessel  in  tissue  (Figure  2  (a)).  A  uniform  compression  in  the 
z  direction  is  applied  and  the  3D  displacement  held  of  the  phantom  is  calculated 
using  ABAQUS  hnite  element  package  (Providence,  RI).  The  Poisson’s  ratio  is 
set  to  u  =  0.49  in  both  phantoms  to  mimic  real  tissue,  which  causes  the  phantom 
to  deform  in  x  &  y  directions  as  a  result  of  the  compression  in  the  z  direction. 

Respectively  5  x  10^  and  1.4  X  10^  scatterers  with  uniform  scattering  strengths 
are  uniformly  distributed  in  the  hrst  and  second  phantom,  ensuring  more  than 
10  scatterers  exist  in  a  resolution  cell.  The  scatterers  are  distributed  in  the  8mm 
diameter  vein  also  (Figure  2  (a)).  To  construct  deformed  ultrasound  images,  the 
displacement  of  all  of  the  scatterers  is  calculated  by  interpolating  the  displace¬ 
ment  of  the  neighboring  nodes  in  the  hnite  element  analysis.  The  parameters 
of  the  probe  are  set  to  mimic  Siemens  5-lOMHz  probes.  The  probe  frequency  is 
7.27MHz,  the  sampling  rate  is  40MHz  and  the  fractional  bandwidth  is  60%. 

The  hrst  phantom  undergoes  uniform  compressions  in  the  z  direction  to 
achieve  strain  levels  of  2%  to  14%  in  2%  intervals.  Ground  truth  integer  dis¬ 
placement  values  are  used  as  the  initial  estimate  for  AM  to  decouple  the  perfor¬ 
mance  of  DP  from  AM.  Accurate  subpixel  displacement  held  is  calculated  with 
AM  and  the  mean  strain  values  are  compared  with  the  ground  truth  (Figure  1 
(a)-(c)).  The  results  are  only  shown  for  2%,  4%,  8%  and  14%  compression  for 
better  visualization.  The  results  with  two  threshold  values  for  IRLS  and  without 
IRLS  demonstrate  that  outlier  rejection  does  not  ahect  the  mean  strain  value, 
while  increasing  the  regularization  weight  increases  underestimation  of  the 
displacement.  The  rate  of  increase  of  the  underestimation  with  increasing  is 
signihcantly  more  with  the  unbiased  regularization  (dashed  line)  as  expected. 

Signihcantly  higher  signal  to  noise  ratio  (SNR)  [2]  values  can  be  achieved 
with  outlier  rejection  (Figure  1  (d)-(f))  without  over-smoothing  the  image  with 
high  Qfa  values.  To  show  the  performance  of  the  overall  method,  the  initial  inte¬ 
ger  displacement  held  is  calculated  with  DP  and  accurate  displacement  held  is 
calculated  with  (Figure  1  (g)-(i)).  The  SNR  values  are  less  than  previous  case 
especially  at  high  strain  values,  where  DP  results  deviates  from  ground  truth. 

The  second  simulation  experiment  is  designed  to  show  the  ehect  of  smooth¬ 
ness  weight  and  IRLS  threshold  on  contrast  to  noise  ratio  (CNR)  [2]  when  the 
correlation  is  lower  in  parts  of  the  image  due  to  huid  motion.  The  phantom 
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(a)  T  =  0.005 


(b)  r  =  0.01 


(c)  no  IRLS 


(d)  T  =  0.005 


(e)  T  =  0.01 


(f)  no  IRLS 


(g)  T  =  0.005  (h)  T  =  0.01  (i)  no  IRLS 

Fig.  1.  Mean  and  SNR  of  the  elastograms  of  the  Field  II  simulated  uniform  phantom 
at  four  different  compression  levels  (shown  in  percentage)  for  three  IRLS  T  values.  The 
solid  and  dashed  lines  correspond  to  biased  and  unbiased  regularizations  respectively. 
(a)-(c)  shows  the  relative  underestimation  of  the  strain,  e  is  the  mean  strain  calculated 
with  the  elastography  method  and  e*  is  the  ground  truth,  (d)-(f)  shows  the  SNR  of 
the  AM.  (g)-(i)  shows  the  SNR  of  the  AM  with  initial  displacements  found  by  DP. 


contains  a  vein  oriented  perpendicular  to  the  image  plane  (Figure  2).  The  initial 
integer  displacement  is  generated  with  DP.  The  background  window  for  CNR 
calculation  is  located  close  to  the  target  window  to  show  how  fast  the  strain  is 
allowed  to  vary,  a  property  related  to  the  spatial  resolution.  The  maximum  CNR 
with  IRLS  is  5.3  generated  at  T  =  0.005  and  eta  =  38,  and  without  IRLS  is  4.8 
at  CKa  =  338.  Such  high  value  makes  the  share  of  the  data  term  in  the  cost 
function  very  small  and  causes  over-smoothing. 

Phantom  Results.  We  perform  freehand  palpation  experiment  on  a  breast 
phantom  to  examine  the  performance  of  the  frame  selection  technique.  50  frames 
of  RF  ultrasound  data  are  acquired  using  a  Siemens  Antares  system  (Issaquah, 
WA).  Our  custom  data  acquisition  program  is  connected  to  the  Axius  Direct 
Research  Interface  to  send  the  command  for  capturing  RF  data.  At  the  same 
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(a)  simulation  phan¬ 
tom 


(b)  finite  element  strain 


(c)  CNR 


Fig.  2.  The  target  (circular)  and  background  (rectangular)  windows  for  CNR  calcula¬ 
tion  of  (c)  are  shown  in  (b) 


Fig.  3.  The  SNR  and  CNR  of  the  phantom  experiment  with  and  without  frame  selec¬ 
tion 


width  (mm) 

(a)  CT  image  (b)  B-mode  image 


width  (mm) 

(c)  strain  image 


Fig.  4.  Patient  experiment  resnlts.  The  arrow  points  to  the  lumpectomy  cavity. 


time,  the  program  collects  tracking  information  from  a  Polaris  tracker  (Waterloo, 
Canada) .  Currently,  the  RF  frames  are  stored  on  the  ultrasound  system  and  are 
processed  offline.  Figure  3  shows  the  SNR  and  CNR  results.  In  automatic  frame 
selection,  Qij  (equation  8)  for  any  two  frames  i,j  in  a  buffer  of  size  15  frames 
is  calculated.  For  the  two  frames  which  give  the  minimum  Q,  the  strain  image 
is  obtained.  The  next  image  is  then  fed  to  the  buffer,  its  first  image  is  removed 
and  the  frame  selection  is  performed  again.  The  automatic  frame  selection  gives 
8  frame  pairs  for  strain  calculation  (as  seen  in  the  figure  by  8  SNR  and  CNR 
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values).  Without  frame  selection,  49  strain  images  are  calcnlated.  The  average 
CNR  and  SNR  values  are  improved  from  4.91  to  7.19  and  from  4.98  to  5.88  with 
frame  selection. 

Patient  Results.  We  have  acquired  freehand  palpation  nltrasound  RF  data 
using  the  Siemens  Antares  system  from  patients  approximately  four  weeks  after 
lumpectomy.  The  ultrasound  probe  is  tracked  with  the  Polaris  tracking  system. 
Optimal  frame  selection  is  performed  to  select  images  for  elastography  nsing  the 
AM  method.  The  strain  image  (Fignre  4)  shows  that  the  AM  method  can  detect 
the  thin  hard  scar  tissne  even  thongh  it  is  close  to  the  cavity  flnids  which  nndergo 
incoherent  motions  and  canse  signal  decorrelation.  Since  the  AM  method  hnds 
the  displacement  of  all  the  samples  on  an  A-line  at  the  same  time,  the  correlated 
data  at  the  top  and  bottom  of  the  cavity  gnide  the  method  to  hnd  the  correct 
displacement  inside  the  cavity  where  the  data  is  decorrelated. 

4  Discussion  and  Conclusion 

We  introduced  a  novel  method  for  calcnlating  a  dense  displacement  map  by 
analytic  minimization  of  a  cost  function.  We  used  the  IRLS  method  from  ro¬ 
bust  statistics  to  make  the  tracking  resistant  to  outliers.  Moreover,  we  exploited 
the  tracking  data  to  optimize  frame  selection.  Through  simulation  studies  using 
Field  II  and  hnite  element  analysis,  we  showed  that  the  proposed  AM  method 
generates  high  qnality  displacement  estimates.  The  elastography  method  works 
in  real-time.  A  comparison  of  the  IRLS  method  with  quality  guided  displacement 
tracking  [14]  which  also  aims  for  robustness  is  a  subject  of  future  work. 

We  chose  the  novel  application  of  the  lumpectomy  cavity  localization  as  the 
hard  scar  tissue  is  relatively  thin  and  demands  a  high  resolution  elastography 
method.  Also,  incoherent  fluid  motions  in  the  cavity  causes  large  decorrelations, 
requiring  a  robust  method.  We  have  an  active  Institutional  Review  Board(IRB) 
protocol  and  have  promising  results  from  9  patients  which  will  be  pnblished  in 
fntnre  work. 
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